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Abstract 

The nonlinear Schrodinger equation is studied for a periodic sequence of delta- 
potentials (a delta-comb) or narrow Gaussian potentials. For the delta-comb the 
time- independent nonlinear Schrodinger equation can be solved analytically in terms of 
Jacobi elliptic functions and thus provides useful insight into the features of nonlinear 
stationary states of periodic potentials. Phenomena well-known from classical chaos 
are found, such as a bifurcation of periodic stationary states and a transition to spatial 
chaos. The relation of new features of nonlinear Bloch bands, such as looped and 
period doubled bands, are analyzed in detail. An analytic expression for the critical 
nonlinearity for the emergence of looped bands is derived. The results for the delta- 
comb are generalized to a more realistic potential consisting of a periodic sequence of 
narrow Gaussian peaks and the dynamical stability of periodic solutions in a Gaussian 
comb is discussed. 

1 Introduction 

In the case of low temperatures, the dynamics of a Bose-Einstein condensate (BEC) can 
be described in a mean-field approach by the nonlinear Schrodinger equation (NLSE) or 
Gross-Pitaevskii equation (see, e.g., [51,52]) 

-^V' + V{x)+g\i.{x,tr^ ^(x,t) = i/,^^ , (1.1) 

where g is the nonlinear interaction strength. Another important application of the nonlin- 
ear Schrodinger equation is the propagation of electromagnetic waves in nonlinear media 
(see, e.g., [23], chapter 8). 
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Due to the possibility to control all experimental parameters accurately over wide ranges 
and monitor the dynamics of a BEC in situ, these have become one of the most prominent 
models for the study of nonlinear dynamical systems. An increasing number of important 
nonlinear phenomena has been demonstrated experimentally in the last decade, such as the 
motion of dark [10] and bright solitons [34] and the self-trapping effect [4]. Furthermore, it 
is possible to reduce the dimensionality of the NLSE in confined geometries (see, e.g., [29] 
and references therein). An excellent review of nonlinear phenomena in BECs can be 
found in the recent article by Carretero-Gonzalez et al. [13]. A particular important 
subdiscipline is the study of the nonlinear dynamics of ultracold atoms and BECs in 
periodic potentials since these systems provide excellent model systems for fundamental 
problems in condensed matter physics (see [45] and [6] for recent reviews). Among other 
phenomena, recent experiments have demonstrated the excistence of gap solitons [24] and 
of looped Bloch bands and the corresponding dynamical instability [25, 33] as well as 
nonlinear de- and rephasing [47, 68] . 

In this paper paper we will present an analytic study of stationary states of the non- 
linear Schrodinger equation, i.e. states of the form ^(x,t) = exp(— i|Ut//i)'(/'(2;), in a one- 
dimensional periodic potential. These states satisfy the time-independent NLSE 



with V{x + d) = V{x) and /i G M. Here we adopt a novel approach to nonlinear stationary 
states by exploiting the analogy to nonlinear dynamical systems in considering the wave 
function ip{x) as a dynamical variable. Although a bit unusual, this approach leads to a 
variety of important results, explaining several phenomena umambigiously in the language 
of nonlinear dynamical systems. For instance, the appearance of period-doubled Bloch 
bands, as described previously in [40], is readily understood as a familiar period-doubling 
bifurcation. Similarly, the bifurcation leading to the appearance of looped Bloch bands is 
easily understood so that it is possible to calculate the critical nonlinearity for this process 
analytically. From a fundamental point of view it becomes clear that Bloch states, i.e. 
solutions which are periodic up to a phase factor e*'"^, are extremely fragile objects. While 
in the linear case g = all states are Bloch states or linear combinations of degenerate 
Bloch states according to Bloch's theorem, these states are of measure zero in the nonlinear 
case g ^ 0- The reason for this difference is that almost all solutions are only quasi-periodic 
for g 0, approaching a true periodicity smoothly in the limit g ^ 0. With increasing 
strength of the nonlinearity, more and more (quasi)-periodic orbits are lost and chaos 
sets in in the wake of a period doubling process. Spatial chaoticity usually implies a 
divergence of the wave function. Thus, well-behaved nonlinear stationary states are very 
rare for certain parameter values, in fact large interaction constants. 

As a convenient model for a space-periodic structure, we consider a particular simple 
model system, the delta-comb potential 



The potential (jl.3p has the great advantage that analytical solutions of the NLSE can be 
constructed in terms of Jacobi elliptic functions [11]. Analytic solutions of the nonlinear 




(1.2) 




(1.3) 
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Schrodinger equation for a non-vanishing potential V{x) are very rare and hence it will 
be of interest to study such a potential in some detail. Its linear counterpart represents a 
basic model for a quantum mechanical Bloch system (see, e.g. [30]). A generalization to 
a more realistic setup consisting of a periodic series of Gaussian potential barriers will be 
discussed at the end of this paper. Narrow Gaussian potential barriers can be realized in 
a good approximation in current atom-chip experiments (see [27] for a recent review). For 
example the resonant transport of BECs through two such peaks is analyzed in [49]. We 
generally restrict ourselves to the case of a repulsive interaction, g > and // > 0. Some 
differences that arise for an attractive interaction will be discussed briefly in Sec. 12. 4[ 
For convenience we set h = M = 1 and fix the period of the potential as d = 2Tr. If 
the normalization of the wave function is arbitrary, one can rescale the wave function 
as V V'/y^; leaving two free parameters (A, /x), otherwise three parameters must be 
considered. 

Up now, quite a number of articles were devoted to the important problem of nonlinear 
stationary states in a periodic potential and we refer to the review articles [13, 45] for 
an overview. There also have been some papers devoted to the NLSE with a delta-comb 
potential, however focusing on completly different aspects. Bound states were studied 
in [60] and the nonlinear transport problem through such a series of barriers was discussed 
in [59]. Bloch bands were studied in [22,38,57,58], however relying mostly on numerical 
approaches. Recently, Paul et al. used a randomized delta-comb potential to study the 
crossover to Anderson localization in [50] in ultracold gases. A related dynamical approach 
to nonlinear eigenstates can be found in [53], however without taking much benefit from 
the theory of nonlinear Bloch bands. It should also be noted that bound states of nonlinear 
Schrodinger equations with a nonlinear periodic microstructure where g = g{x) is space- 
periodic have been explored [26,64]. 

In particular this paper is organized as follows. The stationary solutions of the free 
NLSE (i.e. V{x) = 0) are reviewed in Sec. [2] to lay the foundation for our further work. 
With these results in mind we then construct the solutions for the delta-comb potential in 
terms of Jacobi elliptic function in Sec. 12. 2^ thus reducing the NLSE to a discrete mapping 
for their parameters. In Sec. [3] we focus on periodic solutions and derive conditions for 
their existence. The stability of these solutions is discussed in detail in Sec. 13.21 We then 
apply our results to the important problem of nonlinear Bloch bands in SeclH deriving the 
critical nonlinearity for their emergence analytically. Finally, we generalize our approach to 
a more realistic setup consisting of a periodic series of Gaussian potential barriers in Sec.[5l 
In addition we discuss the temporal stability of nonlinear periodic solutions, showing that 
the stability properties are fundamentally different for different types of nonlinear Bloch 
states. A summary and discussion of our results is presented in Sec. [6] 

2 Solutions of the nonlinear Schrodinger equation 

Before discussing particular solutions, we start with some general remarks on the time- 
independent nonlinear Schrodinger equation (NLSE) 

First, we note that the real-valued time-independent NLSE is mathematically equiva- 
lent to a classical nonlinear oscillator described by the nonlinear Hill equation 



y + /(t)2/ + = 0, f{t + T) = f{t). 



(2.1) 
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Replacing the amplitude y{t) by the wave function and the time t by the position vari- 
able X and identifying (3 = —2g and f{t) = — V(x)) one recovers the time-independent 
NLSE ()1.2|) . For V{x) = one recovers the undamped Ueda oscillator (see, e.g., [36,44,48] 
and references therein). Within the framework of nonlinear eigenstates of BECs, this equa- 
tion has been analyzed in [11]. For a cosine-potential V{x) = Vbcos(x), Eqn. (12. ip is just 
the nonlinear Matthieu equation, which is a popular example in classical nonlinear dy- 
namics (see, e.g., [7,19,43,48]). The NLSE with a delta-comb potential V{x) = X5d{x), 
or analogously the kicked nonlinear Hill equation, studied in the present paper has the 
advantage that it allows analytical solutions. Further exactly solvable examples can be 
found when the potential is given in terms of Jacobi elliptic functions [9,31]. 

In quantum systems, however, we also encounter complex solutions of the NLSE dis- 
cussed in sections H] and [5] in the context of nonlinear Bloch bands and the temporal 
evolution of a wave function, which is governed by the time- dependent NLSE. 

Furthermore, the time-independent NLSE (jl.2p can be written as a Hamiltonian system 
with the conjugate variables (■;/'*, V'O (^jV*'*)) where ip' = dip/dx. Introducing the 
Hamiltonian function 

n = m'-2iV{x)-^iM^-9\ij\\ (2.2) 

Hamilton's equations read 

:f^' = -|^ = 2(l^(x)-^)^ + 25|V'|V = < (2.3) 
dx o^p* 

and analogously for (^/^,^'*). Considering only real- valued solutions, one is left with a 
two-dimensional phase space {ipjip') where the flow is is area-preserving. 
For the special case of a delta-comb, 

V{x) = X62n{x) = X^6{x-27rn), (2.5) 

n 

the solutions are essentially the ones of the free NLSE. Hence we give a brief review 
of the real solutions of the free nonlinear Schrodinger equation in the following section. 
The first derivative of the wave function ^' is discontinuous at x = 27rn (cf. the studies 
[2,20,54,56,67] of a NLSE for a single delta-potential) 

lim^ ('0'(27rn + e) - V'(27rn - e)) = 2AV'(27rn), (2.6) 

whereas the wave function itself is continuous. The discontinuity of the delta-comb poten- 
tial does not affect the area-preserving quality of the flow {ip{xo),jp' (xq)) {ip{x),jp' (x)). 
To clarify this issue, we linearize the flow in the vicinity of a delta peak at x = 27rn : 

/ i;{27Tn + e) \ M 0\f V(2vrn - e) \ 

VV'(2™ + e)y V 2A 1 J \ i;'{2Trn-e) J - ^ ' 



As the determinant of the matrix equals unity, the flow is clearly area-preserving. 
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Figure 1. Solution types of the free nonlinear Schrodinger equation witli a repulsive nonlinearity 
g. Left panel: Dependence of the solution type on the initial condition {^P(xq),'iJj' (xq))- Right 
panel: Examples of the solution types for tp'{xo) — and ip{xo) = 0.8 (sn), 1.2 (dc), 1.8 (nc) and 
M = 5 = 1- 



2.1 Real solutions of the free nonlinear Schrodinger equation 

The free nonlinear Schrodinger equation has well known real solutions in terms of Jacobi 
elliptic functions [11,21] (see, e.g., [1,37] for an introduction). 

It will be instructive to calculate the general solution of the free NLSE explicitly. 
Multiplying Eqn. (jl.2p for V{x) = by -0' and integrating once yields the general solution 
of the free NLSE in the form 

x-xo= / . (2.8) 

The type of solution of Eqn. (|2.8p strongly depends on the initial values '4'{xq) and ip'^xo), 
which also determine the constant of integration E via 

E = iv'(xo)' - |V'(xo)' + f^Hxof. (2.9) 

This dependence on the initial values is illustrated in Fig. [TJ For initial values of ip{xo) 
and iP'{xq) inside the regions marked by sn, dc and nc, the integral in Eqn. (j2.8p can be 
reduced to the canonical form of the Jacobi elliptic functions sn, dc and nc, respectively. 
For, e.g., iP'{xq) = a simple scaling of the wave function 'ip = \/ gj [i^) and the position 
X = ^/Jlx leads directly to the standard forms given in [37], Chapter 3. Examples of these 
solutions types are shown on the right-hand side of Fig. [TJ Note that only the sn-type 
solution is periodic while all other types diverge at a finite value of x. 

The situation is more involved for initial values in the regions marked by G in Fig. [H 
A simple scaling is not sufficient any longer, but a transformation t = {s + q)/{s — q) with 
q = {2E/g)^^^ brings the integral (12. 8p to the standard form of the sc function [37]. 

Note that the sn-type solutions are periodic, while all other solutions (nc, dc and G) 
diverge at a finite value of x. This is easily understood in terms of the classical analog, the 
nonlinear Hill equation (12. ip . This equation describes a classical particle with coordinate 
y and kinetic energy T = y^/2 moving in a potential U{y) = fiy'^ + /3?/^/4 with (3 = —2g, 
that is illustrated in Fig. [2l For g > this potential has a minimum at y = and 
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Figure 2. Classical potential U{y) describing the nonlinear Hill equation rcsp. undamped Ueda 
equation (|2.ip for fj. — 1 and g — I ( — ) resp. g — —1 ( ). 



diverges to — oo for \y\ oo. Thus y will diverge as soon as the particle can leave the 
potential minimum around y = 0. In contrast, the motion is always bounded for g < 0. 
The different regions in Fig. [T] are directly related to the turning points for a system 
with energy E = T + U = + /uy^ + f3y^/A. For g > and energies in the interval 
< E < C/max = /f^^/2(7, we have four real valued turning points with bounded trajectories 
in the sn-region and unbounded ones in the dc-region. For E < 0, we have two real and two 
complex turning points and in the G-region for E > C/max all turning points are complex. 
The curves in Fig. [T]are the boundaries of these regions, i.e. E = C/max and E = 0. 

In the following we focus on the sn-type solutions since they do not diverge. For fixed 
values of fi and g these solutions are explicitly given by 



p], (2.10) 



V " L 

where p {0 < p < 1) denotes the elliptic parameter and L is the period of the wave 
function. The amplitude A and the period L are fixed by 

9(P + 1) (1 

where K{p) denotes the complete elliptic integral of the first kind. However, it should be 
kept in mind that other types of solutions are appropriate for different initial values as 
illustrated in Fig. [H 

2.2 Non-diverging solutions for a delta-comb 

The real solutions of the nonlinear Schrodinger equation for a delta comb are essentially 
the ones of the free equation discussed above. Hence we make the ansatz 



AnSn(^AK{pn)^f^ Pr?j (2.12) 



for X G (27rn, 27r(n + 1)) with An and L„ defined in Eqn. (j2.1ip . However, two conditions 
have to be fulfilled at x = 27rn: The wave function is continuous, whereas its derivative 
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Figure 3. Periodic trajectories are given by solutions of the free nonlinear Schrodinger equation 
connected by vertical jumps caused by the delta potentials. The figure shows the 'configuration 
space' representation (a,b) as well as the phase space representation {ip.ip') (c,d). The 

period-one orbit for A = —0.2 (a,c) bifurcates to a period-two orbit (b,d) if the strength of the 
delta-potentials is changed to A = —0.205. The remaining parameters are g = fi = 1. 



is discontinuous according to Eqn. (j2.6p . This is illustrated in Figl3] (a,b), where the 
discontinuity of the the first derivative at x = 27rn is clearly visible. For the sn-type 
solutions (j2.12p the (dis-) continuity conditions are explicitly given by: 

I. AnSn{Un\Pn) = An+l Sn{Un+l\Pn+l) (2.13) 

II. 2A^„ sn{un\pn) + ^cn(n„|p„)dn(n„|p„) 



L 



n+l 



cn(n„+i|p„+i)dn('u„+i|p„+i) (2.14) 



where the abbreviations ti„ = 4i^(p„)(27r(n -|- 1) -|- Xn)/Ln and Un+i = 4:K{pn+i){2-K{n + 
1) + Xn+i)/Ln+i have been used. If |j4„sn(M„|p„)| < |^n-i-i| the first condition can always 
be fulfilled by an appropriate choice of the "phase shift" Xn^i. Inserting the first condition 
into the second one and using the addition theorems of the Jacobi elliptic functions one 
arrives at 

^"^^ +^^sn^(.„|p„) (2.15) 



{pn+l + 1)2 {Pn + 1)2 Pn + 1 

+ 7 ^" o/n —^Cn{Un\Pn)(in-{Un\Pn)sn{Un\pn)- 

(Pn + ly/^ V2;U 
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These equations define a discrete mapping / : {p„ 



or / : (tpn^tpn) 

{ipn+i,tpn-\-i), respectively, where we introduced the abbreviations tpn = V'(2vrn + 0) and 
V'n — V''(27rn + 0) for convenience. In the following we will mainly consider the mapping 
/ in terms of the physical variables ip and ip' . This is the Poincare section of the area- 
preserving phase space flow generated by the NLSE. However, the mapping / in terms of 
the abstract parameters p„ and Xn is more suitable for actual calculations. 

Note that the discontinuity of the derivative may lead to a divergence of the wave 
function. As stated above, the sn-type solutions are periodic, while all other solutions of 
the free NLSE diverge at a finite value of x. Due to the discontinuity of the derivative, 
the mapping / can map a point {ipn,tpn) inside the sn-region of phase space to a point 
{ipn+i,tpn-\-i) outside the sn-region. This will generally lead to a divergence of ip{x) at a 
finite value of x. We will come back to the question of divergences later in this section. 
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Figure 4. Stroboscopic phase space plots of the mapping / : {tpmi^n) ~^ (V'n+i, V'n+i) for g = 1, 
/.t = 1 and A 0.02, 0.1, 0.5 (from left to right). 



2.3 Dynamics in phase space 

It is instructive to first have a look at the flow in phase space as illustrated in Fig. [3] 
(c,d) for periodic solutions {g = 1, fi = 1) in comparison to the respective wave function 
representation (a,b). Panel (c) show a period-one trajectory following an orbit of the 
free nonlinear system as discussed in Sec. 12.11 up to a point tp « —0.84 where the delta- 
potential causes a vertical jump of ip' according to (j2.6p . If the potential strength A is 
made more negative, the orbit slightly adjusts until it bifurcates into a double-periodic 
orbit as shown in Fig. [3](d) for A = —0.205. Here two kicks connect two different orbits of 
the free nonlinear system. We will come back to the period-doubling bifurcation scenario 
in Sect. 13.21 (compare in particular Fig. [TT]) . Similar phase space plots can be found in 
Ref. [5,32,35]. 

We now turn to the Poincare section of this this flow given by the discrete mapping / or 
/ respectively, because it offers a global view of the dynamics. For the actual computation 
we use the mapping / in terms of pn and x„, in fact the Eqn. ()2.15p . We propagate an 
ensemble of trajectories with randomly chosen initial conditions within the sn-region of 
phase space (cf. Fig. [1]). The dynamics within this region is visualized in Figs. H] and [5] 
for different values of the potential strength A, both for repulsive and attractive delta- 
potentials. The remaining parameters are fixed as fi = 1 and g = 1- First we notice that 
the dynamics is invariant with respect to a global sign change (V'njV'n) ~^ i~'4'n, —i^n)^ 
hence the phase space is inversion symmetric. For A = the phase space is additionally 
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symmetric to the coordinate axes. This symmetry is destroyed with increasing potential 
strength A. 
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Figure 5. As figure[4l however for A = —0.02, —0.1, —0.5 (from left to right). 



The dynamics is quasiperiodic for small amplitudes, i.e. in the vicinity of the trivial 
solution {^p{x),^p' (x)) = (0,0). It is chaotic for large amplitudes at the edge of the sn- 
region, corresponding to a stronger nonlinear interaction. The chaotic part of phase space 
becomes larger with increasing potential strength |A|. The qualitative difference between 
stability and chaos is illustrated in Fig. [6l the 'dynamics' of the wave function ^(x) in 
configuration space and in phase space for n = g = 1 and A = —0.2 and A = —0.22, 
respectively. In both cases, the initial state V'(O), '(/''(O) was chosen as a periodic solution 
plus a small random perturbation (standard deviation a = 10~^). This periodic solution is 
elliptically stable for A = —0.2 so that the perturbation does not grow with x. It becomes 
hyperbolically unstable for A = —0.22 so that the perturbation increases exponentially. 
Trajectories with different perturbations spread rapidly, as shown in Fig. [6] (c) and (d). 
The different forms of fixed points and their stability in dependence on the parameters 
will be discussed in detail in the next section. 

As illustrated in Fig. [6] (c), a trajectory in the region of spatial chaos will generally 
diverge at some finite value of x. Such a trajectory is usually mapped close to the edge of 
the sn-region in phase space, where non-diverging solutions exist (cf. Fig. [T]). Due to the 
discontinuity of the derivative ip'{x) the wavefunction leaves the sn-region. As mentioned 
above, all non-sn-type solutions of the NLSE diverge at a finite value of x. Thus also 
trajectories started in the chaotic part of the sn-region usually diverge at a finite value of 
x. An example of such a diverging wavefunction is shown in Fig. [6] (c,d) for ^ = g = 1 
and A = —0.22 (cf. Fig. [5|). The initial value {ip{0),'ip' (0)) was chosen as an unstable fixed 
point of / inside the chaotic region plus a small random perturbation. 

The stability or divergence of a wavefunction depends sensitively on the initial values 
{il>{0),'ip' {0)). This is illustrated in Fig. [71 A grey-scale plot shows the point of diver- 
gence (the number of spatial periods until a trajectory diverges) in dependence of the 
initial condition {iPq,iP'q) for A = 0.5. The left figure shows the position of divergence 
for one quadrant of phase space, ipo,ipQ > 0. This figure should be compared with the 
Poincare section in Fig. [H One clearly recognizes that trajectories with small amplitudes 
are quasiperiodic and thus stable. But one can also find initial values that do not lead to 
divergences for larger amplitudes. These values form an approximately self-similar set in 
phase space. This is illustrated in Fig. [7| in the lower plots, where magnifications of the 
upper plot are shown. 
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Figure 6. (In)stability of the wave function t/j{x) for ji = g = \ and A = —0.2 (a,b) and A = —0.22 
(c,d). Tlie plots shows tlie 'dynamics' of the wave function ipi^) in configuration space (a,c) and in 
phase space (b,d). In all cases, the initial state V(0),V''(0) was chosen as the symmetric periodic 
solution plus a small random perturbation (standard deviation a = 10~^). The periodic solution 
is elliptically for stable A = —0.2 (a.b) so that the error docs not grow and the wave function i!{x) 
remains close to it for all x. The periodic solution becomes hyperbolically unstable for A = —0.22 so 
that the error grows exponentially and the the wave function eventually diverges. Ten trajectories 
with different random perturbations are plotted in the figure, showing the rapid spreading of these 
trajectories. 

2.4 Attractive nonlinearity 

For the sake of completeness we also briefly discuss the case of an attractive nonlinearity 
g <Q without going into details. 

Real solutions of the nonlinear Schrodinger equation for a delta comb can be constructed 
in terms of the Jacobi elliptic function cn and dn. The cn-type solutions are given by 

Pn) (2.16) 

for X G (27r?7., 27r(n + 1)) with 



V'(a;) = AnCn\^K{pn) 



X + Xr. 



(2.17) 
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Figure 7. Number of spatial periods, at which a trajectory diverges in dependence of the initial 
condition (ipoii^'o) ^^r A = 0.5. The lower panels show subsequent magnifications. 



The dn-type solutions read 



Lr 



Pn 



for X £ {2-1171, 27r(n + 1)) with 



fi(2-p) 



and 



ijp - 2)K{pf 



(2.18) 



(2.19) 



Note that cn-type solutions exist for /i < 0, leading to p G (0.5, 1], as well as for /U > 0, 
leading to p G [0,0.5), whereas the dn-type solutions exist only for /i < 0. 

Again, two conditions have to be fulfilled at x = 2TTn: The wavefunction is contin- 
uous, whereas its derivative is discontinuous according to Eqn. (j2.6p . As above we can 
consider the dynamics stroboscopically at and thus arrive at the mapping / : ~^ 
(V'n+i) V'n+i)- Examples of the dynamics in phase space in a delta comb of strength 
A = -0.1 are plotted in Fig. El 

An important difference to the case of a repulsive interaction is that one can find a 
solution of the free NLSE in terms of the non-diverging Jacobi elliptic functions cn and dn 
for all initial values of a wave function tp{xo),^p' (xq). This should be compared to the case 
of a repulsive interaction, in particular to Fig. [TJ Hence the solutions for a delta-comb do 
not diverge at a finite value of x, regardless of the starting point (V'Oi V'o)- The dynamics in 
phase space around {^n,ipn) = (0,0) is illustrated in Fig. [8]for A = —0.1 and two different 
values of 
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Figure 8. Stroboscopic phase space plots of the mapping / : (?/;„, '0^^) {ijjn^i^ilj'^-^j^i) for an 
attractive nonhnearity g — —1 and fi — +1 (left panel) and fi — —2 (right panel). 

The solutions are always given by the Jacobi elliptic function cn for fi > 0. For /x = 1 
and A = —0.1, the dynamics is still completely regular and one we can identify fixed points 
again (cf. the left panel of Fig. [8|). The dynamics becomes partly chaotic for higher values 
of A, but still a wave function cannot diverge at a finite value of x. However, a totally 
different mechanism of divergence exists in this case. The wave function tp{x) diverges 
for x ^ oo, if its phase is such that the changes of the amplitude due to the delta kicks 
accumulate. Revisiting Eqn. (j2.17p one recognizes that this divergence of the amplitude 
An implies pn — > 0.5. Indeed one finds numerically that most trajectories in the chaotic 
region of phase space show such a behavior. 

For fj, < the solution is described either by the function cn or dn, depending on the 
initial value ('i/'OjV'o)- '^^^ dynamics in phase space is illustrated in the right panel of 
Fig. [8] for A = —0.1. Solutions of dn-type are found around the elliptic fixed points at 
(ipm^'n) = (il) 0). One can identify regular islands in a chaotic sea, surrounded by quasi- 
periodic orbits for large amplitudes. Again the cn-type solutions can diverge for x ^ oo 
if Pn 0.5. 



2.5 Complex solutions 

Complex solutions of the free NLSE can be found by decomposing the wave function into 
amplitude and phase 



with real-valued functions (pix) and S{x) > 0. The phase is then given by 



<P'ix) 



a 



S{x) 



(2.20) 



(2.21) 



with a real integration constant a. In the limit a — > one recovers the real- valued sn- 
and cn-type solutions described above, which change sign at the zeros of the density. This 
limit is discussed in detail in [55] where interaction-induced transitions from scattering 
states to bound states are analyzed for a finite square well potential. 
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We focus on the special solution [11,55] 

Six)=Bn-^dn^(^Qnix + Xn)Pn^ (2.22) 

of the free NLSE with Qn = ^K{pn) / Ln. Inserting this ansatz into the free NLSE yields 
the following conditions for the remaining constants 

2^ = "igBn - {2-pn)Qn 

QtiPn - l)^n + gal + 2glBn - 2gfiBl = 0. 

For the delta-comb potential, one can make the ansatz (12.221) separately for each interval 
X £ (27rn, 27r(n + 1)). The derivative ip' of the the wave function is again discontinuous at 
X = 2irn according to Eqn. (12.61) . In terms of the density S{x) this yields 

lini^ (S'(27rn + e) - S'{2TTn - e)) = 4A5(2^n). (2.23) 

Complex periodic solutions of this type will be discussed in Sec. 13.31 These solutions are 
of great importance, since they occur for nonlinear Bloch bands (cf. Sec. S]). 



3 Periodic solutions 

Real periodic solutions can be viewed as fixed points of the mapping / resp. / defined in 
Sec. 12.21 In the following we calculate these solutions explicitly and discuss their spatial 
stability. 

3.1 Symmetric and antisymmetric real periodic solutions 

One kind of periodic solutions is found for ip{2TTn) = 0. They will be named antisymmetric 
periodic solutions in the following since they fulfill tp{2'im — x) = —'ip{27rn + x). In this 
case the wave function and its derivative are continuous everywhere. The period length of 
the Jacobi elliptic functions L must fulfill 

L = 47r/m (3.1) 

for m G N. The antisymmetric solutions are 27r-periodic only for even m, since their 
fundamental period is Air/m. Note that m is the number of nodes of the wave function in 
[0, 27r). In the following we will mostly consider the "primary resonance" m = 2 with the 
fundamental period 27r. However, we will show in Sec. 13.21 that this solution bifurcates 
when A is varied and new period doubled solutions emerge. 
In terms of the elliptic parameter p Eqn. (|3.ip reads 



{p + l)K{pf = ^. (3.2) 

The left side of this equation is bounded from below by {p + 1) K{p)'^ > i^(0)^ = 7r^/4. 
Thus antisymmetric periodic solutions only exist if /i > Comparing with the case of 
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(b) 2 
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1 l-S 

x/271 





Figure 9. Periodic solutions i^ix) of the NLSE for a delta-comb potential of strength A = 1 for 
/Lt = 3 and 5 = 1. The red vertical lines indicate the positions of the delta potentials. Antisymmetric 
solutions (a,c) have a node at the positions of the delta potential, iJ;{2TTn) — 0, so that they 
experience no discontinuity of the first derivative and are thus independent of A. Symmetric 
solutions (b,d) are non-zero at a; = 27rn and symmetric with respect to a reflexion at this point. 
Furthermore the solutions are characterized by the number of nodes in the interval [0, 27r) - the 
examples have m = 2 (a,b) and m = 4 (c,d) nodes, respectively. 



a harmonic potential, i.e. the nonlinear Mathieu equation [43,48], one recognizes that an 
antisymmetric periodic solution comes into being just at the 2 : m resonance 2/i = m? /A. 

The shift Xn is given by = resp. = vr. However, these solutions differ only by 
an overall sign. Note that the solution of Eqn. (j3.2p does not depend on A. As an example 
the antisymmetric periodic solutions with m = 2 and m = 4 re shown in Fig. [9] (a) and 
(c), respectively, for |U = 3 and g = 1. 

Furthermore, 27r-periodic solutions can be found which are symmetric around the po- 
sitions of the delta potentials at x = 27rn, n G Z. This symmetry and the periodicity of 
the wave function imply that the solution is also symmetric around x = (2n -|- l)7r, as one 
can easily see: 

'4}{{2n + l)7r -x)= '4j{-{2n + 1)tt + x) = V'((2n + 1)tt + x). (3.3) 
Hence the the wave function assumes a maximum or minimum at x = (2n -|- l)7r, in the 
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middle between two delta potentials. Without loss of generality we conclude that 

x„ = -(2n + l)7r + L/4 or x„ = (2n - l)7r + 3L/4. (3.4) 

These solutions differ only by an overall sign, hence only one of them must be considered. 
Again condition (j2.6p must be fulfilled at the positions of the delta-comb. Because of the 
periodicity it is, however, sufficient to evaluate this equation at x = only. A further 
simplification arises from the symmetry of the wave function, yielding the condition 

4:K(p) 

Xsn{u\p) = — - — cn{u\p) dn{u\p) (3.5) 

with u = 4:K{p)xo/L = K{p){l — Air/L). Eqn. (j3.5p has a solution for every n > 0. For 
higher values of fi it may have several solutions, differing in the period L and correspond- 
ingly in the number of nodes of the wave function. Again, these solutions will be labeled 
by the number m of nodes of the wave function in the interval [0, 27r). In the following we 
consider the solution with m = 2 nodes. Exemplarily the solutions with m = 2 and m = 4 
are shown in Fig. [U] (b) and (d), respectively, for /i = 3, 5 = 1 and A = 1. 

One can easily see that these two solution classes are the only possible 27r-periodic 
sn-type wave functions. Due to the periodicity the elliptic parameter pn is the same for all 
n, i.e. pn = P- As the elliptic parameter p is fixed, the derivative ip'{xo) at a point xq is 
determined by the wave function 'ip{xo) up to a sign. Considering xq = 2Tm this leaves the 
possibilities ip'{2Tm + e) = 'ip'{27rn — e), for which one recovers the antisymmetric periodic 
solutions, whereas the other possibility ijj'{2TTn + e) = — V''(27rn — e) leads to the symmetric 
periodic solutions. 

3.2 Spatial Stability of the real periodic solutions 

As already stated, real periodic solutions of the NLSE in a delta-comb are the fixed points 
of the mapping / resp. /, defined in Sec. 12.21 Similarly, we can find solutions with a 
periodicity of 27rr as fixed points of the mapping In the following we discuss the 
characteristics of these fixed points, especially their spatial stability, in dependence of the 
control-parameter A. 

First we discuss the antisymmetric periodic solutions with L = 27r/m for a repulsive 
nonlinearity g = 1 and fi = 1. As an example, we consider the solutions with m = 2, 
which have nodes at x = vrn, n G Z, i.e. at the positions of the delta-potentials x = 2Tm 
and additionally in the middle between them. These periodic solutions are found for 
p PS 0.4719 and x„ = 0, resp. Xn = vr, differing only by an overall sign. 

To determine the stability of a fixed point numerically we iterate the mapping /, where 
the initial value ipo,xo) is chosen as the fixed point plus a small random perturbation. 
Furthermore we calculate the stability (or monodromy) matrix of the mapping / 



(3.6) 



at the fixed point (tp^,tp'^). As the mapping / is area-preserving, the product of the 
eigenvalues exp(ib7) of M is unity. The mapping is unstable if the stability exponent 7 is 
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Figure 10. Phase space around the antisymmetric fixed point of / (marked by a +) before the 
bifurcation (A = 0.14, left) and after the bifurcation (A — 0.16, right). 




Figure 11. Period doubhng of the antisymmetric fixed point with m = 2 (left panel) and of the 
symmetric fixed point (right panel) for /x = = 1. The Fixed points of /, and are plotted 
as a function of the control-parameter A. Stable fixed points are indicated by a solid line, unstable 
fixed points by a dashed line. 



real. Otherwise it is purely imaginary and the mapping is elliptically stable. Then |7| is 
called the stability angle. 

Exemplarily we consider the antisymmetric fixed point with m = 2, which is found 
to be elliptically stable if A > and A is not too large. The stability exponents 7± are 
purely imaginary and one finds quasiperiodic orbits around the fixed point. This is shown 
in the left panel of Fig. [10] for A = 0.14. The fixed point becomes unstable when the 
control-parameter A is increased above a critical value Ac ~ 0.156. This is illustrated in 
the right panel of Fig. [10] for A = 0.16. 

In fact, a bifurcation occurs at the critical value Ac and two new fixed points of 
emerge. This is illustrated in Fig. \\\\ where t/'n = 'i/'(27rn) is plotted for the fixed points of 
Z'', r = 1, 2, 4 as a function of A. One clearly observes a pitchfork bifurcation scenario. The 
new fixed point of no longer fulfills ■0(27rn) = 0, but instead '^(2'Kn) = — ?/)(27r(n + 1)). 
A second bifurcation occurs at A^ ~ 0.178, where the fixed points of become unstable 
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Figure 12. Stability exponent Re(7) for the antisymmetric fixed point witli m = 2 (two nodes 
in [0, 27r), upper panel) and the symmetric fixed point (lower panel) as a function of /i and A in a 
greyscale plot. Wave functions corresponding to these fixed points are illustrated in Fig.[9l Stable 
regions are colored in white (Re(7) = 0). 

and four new elliptically stable fixed points of emerge. These fixed points again become 
unstable at A^' « 0.182. 

One should keep in mind that the NLSE represents a Hamiltonian system and that the 
period doubling scenario of Hamiltonian systems (cf. [3,42]) is different from the familiar 
Feigenbaum scenario for dissipative systems. For instance, the stable fixed points are 
always elliptically (resp. neutrally) stable and not asymptotically stable. 

For A < the antisymmetric fixed points are unstable. However, large regular regions 
exist in phase space around the trivial solution ^(x) = and the symmetric fixed points 
as long as |A| is not too large (cf. Fig. [5]). A "trajectory" started in the vicinity of an 
antisymmetric fixed point is unstable but it may get trapped in the chaotic sea between 
the large regular regions and consequently show quasi-regular dynamics over many periods. 
For stronger potentials, i.e. larger values of |A|, the regular regions shrink and cannot trap 
such a trajectory any longer so that it will in general diverge. 

A similar stability behavior is found for the symmetric fixed point. However, this solu- 
tion is always unstable for a repulsive delta-comb, A > 0, and shows a bifurcation scenario 
for an attractive delta-comb, A < 0. For small repulsive potentials, a similar trapping 
scenario occurs as described above for the antisymmetric fixed points in an attractive 
delta-comb. Again this trapping is lost for stronger potentials, since the regular regions 
of phase space shrink and the trajectories finally diverge in general. 

The symmetric fixed points are elliptically stable for an attractive delta-comb as long as 
A > Ac . A bifurcation occurs if A is decreased below the critical value Ac ~ —0.203. The 
47r-periodic states are again elliptically stable up to the next bifurcation. This bifurcation 
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Figure 13. Stability map of the periodic states for a delta-comb with g = 1. The periodic 
solutions are stable in the shaded areas of the parameter plane (antisymmetric states for A > 
and symmetric states for A < 0). The different solutions are labeled by the number of nodes m in 
the interval [0, 27r). 

route is illustrated in Fig. [11] for the symmetric periodic solution with m = 2 nodes in the 
interval |0, 27r). 

Note that also the trivial solution ^(x) = of the NLSE shows a bifurcation when A is 
decreased. The trivial solution is stable for A > A^^^^ ~ —0.391, where two new elliptically 
stable period doubled fixed points emerge. Quasiperiodic orbits around the emerging fixed 
points are illustrated in the lower panel of Fig. [5] for A = —0.5. 

The stability of the antisymmetric and the symmetric fixed point with m = 2 is sum- 
marized in Fig. [12l where the stability exponent Re(7) is plotted as a function of fj, and 
A. Fig. [13] shows a stability map of several fixed points with different m. 

One can easily understand the qualitative differences of the stability between A < 
and A > by considering the full phase-space dynamics of the classical analogon (|2.ip . An 
antisymmetric periodic solution is not affected by the delta-comb because of '0(27rn) = 0. 
However, a trajectory with a slightly larger classical energy moves ahead of the periodic 
solution, hence it experiences a kick. This kick will increase the energy of the classical 
oscillator if A > 0. The effect of the kicks accumulate and consequently the trajectory 
becomes unstable. If A < the kick will lower the energy and stabilize the trajectory. The 
contrary effect occurs for the symmetric periodic solution. However, the situation is a bit 
more involved here since the periodic solution itself experiences the delta-kicks. 

3.3 Complex periodic solution 

Complex periodic solutions can be found using the ansatz (|2.22p . By similar arguments 
as above one finds that the density S{x) = \ip{x)\'^ must be symmetric around x = 27m. 
Hence it must assume a maximum or minimum at x = (2n -|- 1)tt, in the middle between 
two delta-peaks. Therefore the "phase shift" x^ is given by 




— (2n-|-l)7r or x„ = 



-(2n + l)7r + L/4. 



(3.7) 
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Using the symmetry of the density, the continuity of Eqn. ()2.23p yields the condition 

— — sn(n|p)cn(n|p)dn(ti|p) = B dn^fnlp) (3-8) 

g 

with u = 4:K{p)xn/L for a periodic complex solution. An example of such a periodic 
solution is shown in Fig. [15] (dashed line). 

For a — > the complex solutions tend to the real-valued sn-type solutions which give 
rise to antisymmetric periodic solutions as discussed above. In fact, this is exactly what 
happens at the occurrence of loop structures in nonlinear Bloch bands at /« = (cf. Sec. [4]). 



4 Nonlinear Bloch bands 

Recently, nonlinear Bloch states and Bloch bands have attracted considerable attention. 
New features such as looped Bloch bands [65,69,71] and period doubled Bloch states [40] 
were found. Nonlinear Bloch bands for the delta-comb were first calculated by Seaman et 
al. [57]. In this section we want to link these results to the properties of periodic solutions 
discussed in Sec. [3l 

The appearance of looped bands, or more generally the bifurcation of eigenvalues in 
nonlinear systems, is intimately related to branch point singularities, i.e. exceptional 
points (see [14] and references therein). Those exceptional points are eigenvalue and 
eigenvector degeneracies, which for linear systems can only appear in a non-hermitian 
description. As a boundary of regions of purely real eigenvalues they appear in so called 
PT-symmetric models. In this context, bifurcation phenomena of Bloch bands are also 
discussed in linear, however non-hermitian periodic potentials with PT-symmetry [17, 
18,41]. A full understanding of the interrelation of these phenomena has not yet been 
achieved and may be provided by an analysis of the full nonlinear and non-hermitian 
periodic potential in an extension of the analysis of the corresponding two-state case (see, 
e.g., [28,66]). This is, however, only in the beginning (see [46] for a very recent study). 

Looped Bloch bands arise in nonlinear systems when a space-periodic stationary so- 
lution of the NLSE (a Bloch state with quasimomentum k = 0) undergoes a bifurcation 
associated with the emergence of a dynamical instability. Before the bifurcation, this 
Bloch state is neutrally stable with respect to small perturbations in the sense that the 
time-dependent wave function 'ip{x,t) will remain close to the initial state ip{x,t = 0). 
In the bifurcation, two novel neutrally stable Bloch states emerge, while the old state 
becomes hyperbolically unstable in the the sense that the time-dependent wave function 
ip{x, t) rapidly spreads from the the initial state ^(x, t = 0) leading to full spatio-temporal 
chaos, an example of which is shown in Fig. [22j Such a dynamical instability leads to a di- 
vergence of the Bose-Einstein condensate mode so that a simple mean-field approximation 
is no longer valid. Nevertheless, the time-dependent NLSE is still meaningful since it accu- 
rately predicts the occurence of the instability and the rate of the depletion [15,16,62,63]. 

On the contrary, period-doubled Bloch bands come into existence when a usual 2ti- 
periodic solution ijj{x) of the NLSE undergoes a period-doubling bifurcation to form a 
47r-periodic solution. Consequently, the 27r-periodic wave function is spatially stable up 
to the bifurcation where it becomes spatially unstable as illustrated in Fig. [6l 
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4.1 Linear and nonlinear Bloch states 

Bloch states are nonlinear eigenstates of the form 

i,^{x) = e''^"tx«(x), (4.1) 

where u^{x) is 27r-periodic and k denotes the quasimomentum. The periodic functions 
Uk{x) fulfill the differential equation 

1 / d \ ^ 

~2(dx~^^'*) '^'^^^) + '^^^)'^'^^^) + 9\uk{x)\'^u^{x) = ^x{K)ut,{x) (4.2) 

The eigenenergies ^{n) form the nonlinear Bloch bands. In the preceding sections we 
analyzed the solutions of the NLSE in dependence of the parameter ^, whereas the nor- 
malization of the wave function was arbitrary. For example, periodic states, which are 
Bloch states for k = 0, exist for any value of /i > (cf. Sec. 13. ip . Considering /i as a 
parameter is of course not appropriate here, since this would lead to a continuous and 
multiply degenerate spectrum of //(k) for every k, whereas the corresponding Bloch states 
differ by their normalization. Consequently we consider 27rA^-periodic Bloch states here, 
whose normalization is fixed as 



1 



2'kN 



2ttN 

|u«(x)|2dx = l (4.3) 



throughout this section, leading to a discrete Bloch spectrum. The Bloch states and the 
bands are calculated as described in [71]. 

The set of Bloch states is characterized by one continuous parameter, the quasimo- 
mentum K, and an integer number counting the bands. General non-diverging stationary 
states are given by four continuous parameters, the real and imaginary parts of the initial 
value (V'O) V'o)- Thus the Bloch states are of measure zero within the the set of all possible 
non-diverging stationary states. The nonlinear stationary states in a periodic potential 
generally do not fulfill Bloch's theorem. 



4.2 Looped Bloch bands 

The nonlinear Bloch bands for a delta-comb potential V{x) = \62-k{x) are shown in Fig. [TH 
for (7 = l/vr and A = 0.5. The left-hand side shows the lowest bands in comparison to the 
linear case g = 0, whereas the right hand side illustrates the emergence of loop structures 
at the upper edge of the lowest band at k = 0.5. Similarly loops emerge at the edge of the 
first excited band at «; = 0. The chemical potential is larger than in the linear case by an 
amount of approximately 

\u^{x)\^dx (4.4) 

because of the repulsive mean-field potential. Considering the energy per particle instead 
of the chemical potential one finds swallow's tail structures instead of the loops [57,65]. 

The Bloch states with k = are strictly periodic and can be chosen real up to the 
emergence of loops at a critical nonlinearity ^loop- These states can be identified with the 
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Figure 14. Nonlinear Bloch bands and emergence of loop structures for a delta-comb potential. 
Left: The lowest nonlinear Bloch bands for g — l/n (solid lines) in comparison to the linear case 
g = (dashed lines). Right: The lowest Bloch band for 2TTg = 0, 0.5, 1, 1.5, 2 (from bottom to top) 
and A — 0.5. 



symmetric resp. antisymmetric periodic states described in Sec. 13. li It is found that the 
Bloch state in the lowest band is a symmetric periodic solution with no node, while the 
state in the first excited band is an antisymmetric periodic solution with m = 2 nodes in 
the interval [0, 27r). 

Loops emerge for g > gioop; hence additional Bloch states with k = come into being. 
A state corresponding to the self-crossing of a looped band at k = cannot be chosen 
purely real, whereas the state at the top of the loop can. It is found that the state at the 
top of the loop in the first excited band still corresponds to a real antisymmetric periodic 
solution as for g < ^loop- The Bloch state at the self-crossing corresponds to a complex 
periodic solution as introduced in Sec. 13. 31 The Bloch states in the two lowest bands at 
K = are shown in Fig. [15] for g = l/ir and A = 0.5. The squared modulus of a symmetric 
periodic solution (ground band) and the antisymmetric periodic solution with m = 2 (first 
excited band) are plotted as a bold (symmetric) resp. solid line (antisymmetric). The 
complex solutions that degenerate at the self-crossing at k = transform into each other 
by a sign change of k and consequently by a sign change of x (cf. Eqn. ()4.2|) ). The squared 
modulus of these two solutions is the same and hence symmetric around x = (dashed 
line in Fig. 1151) . whereas the phases are antisymmetric: ipi—x) = —(j){x). 

Using the explicit form of the real and complex periodic solutions one can derive the 
critical nonlinearity g\oop for the emergence of loops in the first excited Bloch band. To 
this end we reconsider condition (jS.Sp for the complex periodic solutions. For g g\oop 
the complex periodic solution and the real antisymmetric solution merge. In this limit 
the wave function tends to zero at the positions of the delta-peaks, ■0(27rn) — > 0, and its 
derivative becomes continuous. Hence both sides of Eqn. (13. Sp tend to zero. 

Consider a complex periodic solution slightly above the critical point. The "wave 
number" q is written as g = Q\oop + ^Q, where 0oop = mK{p)/TT is the wave number at 
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Figure 15. Left: Squared modulus of the wave function of nonlinear Bloch states for k = 0, 
g = \/tt and A = 0.5 (bold line: lowest band, dashed line: state at the self-crossing of the loop 
in the first excited band, thin solid line: state at the top of the loop). Right: Dependence of the 
critical nonlinearity gioop for the emergence of loops in dependence of the potential strength A 
(solid line) and the approximation (|4.9|) (dashed line). 



the critical point g = 5ioop- Expanding Eqn. (|3.8p around 0oop with fixed B then yields 



loop 



9 



t:5q + 0{5q^) = 2XB 



2Xq 



loop 



oop 



9 



9 



Sg + 0{Sg'^ 



(4.5) 



Comparing the coefficients of the linear term yields an equation for the elliptic parameter 
Pioop at the critical point 

2ttX 



PloopKipioopY 



The normalization of the wave function at the critical point yields 



27r 



2n 



sn^ 4K(p) 



L 



P 



2m^ 
■Kg 



Kip) [Kip) - Eip)] 



(4.6) 



(4.7) 



where Kip) and Eip) denote the complete elliptic integrals of the first and second kind, 
respectively. Hence the critical nonlinearity is given by 

,2 



ffloop 



m 



^(Ploop) [-^^(Ploop) - Eipioop)] 



(4.8) 



where pioop is determined by Eqn. (14. Gh . One can derive an approximate explicit expression 
for the critical nonlinearity parameter ^loop by expanding the left-hand side of Eqn. (j4.6p 
up to second order in pioop and substituting the result into Eqn. ()4.8p . This yields in 
second order in pioop 



S'loop 



A 



+ 



A2 



The dependence of 5ioop on A is shown in the lower panel of Fig. [15] for m 
corresponds to the first excited Bloch band. 



(4.9) 
2, which 
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Figure 16. Left panel: Bifurcation of the 27r-periodic state corresponding to the first excited 

Bloch band by a variation of g. Shown is the chemical potential /z of the 27r-periodic state ( ) 

and the emerging 47r-periodic state ( — ). Right: Wave function of the 27r-periodic ( ) and the 

47r-periodic ( — ) state for g = 1/8. 



4.3 Period doubled Bloch bands 

Let us finally discuss another new feature of nonlinear Bloch bands - the existence of 
period doubled Bloch bands [40]. These Bloch bands come into being if a periodic state 
- a Bloch state with k = - becomes spatially unstable in a period doubling bifurcation, 
giving birth to a 47r-periodic state. This state is now again embedded into a (period- 
doubled) Bloch band. 

The period doubling bifurcation has already been introduced in Sec. 13.21 for the case 
of a variation of the potential strength A. As shown on the left-hand side of Fig. \TT\ the 
antisymmetric 27r-periodic state bifurcates when A is increased above a critical value Apd- 
(We now use the index 'pd' for period doubling to distinguish it from the critical value for 
the emergence of looped Bloch bands). 

The period-doubling bifurcation can also take place if the nonlinearity g is increased 
while A and the normalization (14. 3p remain fixed. This is illustrated on the left-hand 
side of Fig. [TBI where we have plotted the chemical potential ;U as a function of g for the 
27r-periodic state (green) and the 47r-periodic state (red). The wave functions of these two 
states are compared on the right-hand side of Fig. [T6j for g = 1/8. 

The asymmetric 27r-periodic state is just the Bloch state with k = in the first excited 
band which is shown in Fig. [T7] (green line). Similarly, the 47r-periodic state is embedded 
in a 'period-doubled' Bloch band which is also shown in Fig. [T7] as a red line. 

The connection to the period-doubling bifurcation of an antisymmetric periodic state 
has several implications for the existence and stability of nonlinear Bloch states: Firstly, 
period-doubled Bloch bands obviously do not exist for g = 0, they come into being when 
g is increased over the critical value i^pd for a period doubling bifurcation. Secondly, it 
was already argued in Sec. 13.21 that the corresponding 27r-periodic state becomes spatially 
unstable in this bifurcation. The lowest symmetric periodic state corresponding to the 
ground band is always unstable for A > 0, while the anti-symmetric state corresponding 
to the first excited band is is always unstable for A < (cf . the stability map in Fig. [T3]) . 

This has several implications for the existence and stability of nonlinear Bloch states. 
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Figure 17. A 47r-periodic Bloch band ( — ) around k = and the corresponding corresponding 
first excited 27r-periodic Bloch band ( ) for A — 0.7 and g— 1/8. 

The 47r-periodic states exist only after a bifurcation, i.e. for g > g^^. They are stable until 
the next bifurcation. As argued in Sec. 13.21 the corresponding 27r periodic state becomes 
spatially unstable at this bifurcation. 

5 The Gaussian comb 

Finally we present some numerical results for a more realistic potential, consisting of a 
periodic structure of narrow Gaussian peaks, 



with Ax = 0.1. The area of each Gaussian peak is normalized to unity. The paramters 
are such that Ax <C ^ <C where ^ = 1/ ^figp is the healing length and p is the averaged 
condensate density. This means that changes in the condensate density are 'healed' on 
a length scale which is much longer than the width of the Gaussian potential barriers, 
but smaller than the separation of the barriers. This gives a rough argument why the 
Gaussian comb is indeed well approximated by a delta-comb discussed above. Such smooth 
potentials can be implemented, at least in principle, in current atom-chip experiments 
(see [27] for a recent review, see also the discussion in [49]). 

5.1 Solutions of the time-independent NLSE 

We consider the solutions of the time-independent nonlinear Schrodinger equation in a 
spatial Poincare section in phase space as in Sec. [21 The Figs. [TSl and [T9l show {ipn,tp'n) = 
(-0(27™), ■0'(27rn)) obtained by a numerical integration for different values of the potential 
strength A. The dynamics in phase space strongly resembles the one for the delta-comb 
illustrated in the Figs. S] and [5j 

Again we find trivial and symmetric periodic solutions, whose spatial stability proper- 
ties are very similar to the case of a delta-comb. A pitchfork bifurcation of these periodic 
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Figure 18. Poincare section in phase space {ipn, ij'n) for the NLSE in a Gaussian comb with g = 1, 
^ = 1 and A = 0.02, 0.1, 0.5 (from left to right). 
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Figure 19. As figure[T8l however for A = —0.02, —0.1, —0.5 (from left to right). 

solutions occurs when a control-parameter is varied. The bifurcation of the trivial peri- 
odic state for A > is illustrated in the upper panel of Fig. [20l The lower panel shows a 
47r-periodic wavefunction that emerges in the bifurcation. 



5.2 Dynamical (in) stability 

Up to now only the spatial stability of periodic solutions has been discussed, not consider- 
ing the temporal behavior, induced by the time-dependent nonlinear Schrodinger equation 



1^ 



+ Vix)+g\i;{x,t)\^]^{x,t) 



dip 



(5.2) 



In the following we analyze the temporal stability of the periodic solutions in the Gaussian 
comb discussed above. The time-dependent NLSE can lead to classical chaotic dynamics 
[61] and dynamical instability, which usually leads to a depletion of the Bose-condensed 
phase [15]. Related work on the dynamical stability of the NLSE in a periodic potentials 
can be found in [8,9, 12]. 

The energetical stability of a Bloch state Mk(x) can be calculated from perturbation 
theory [70,71]. We make the ansatz 



iIj{x) 



[ut,{x) + 5uk(x)] 



(5.3) 



with a small perturbation 5uti{x), that can be decomposed into different modes e^"'^. A 
Bloch state Uk{x) is energetically stable and hence represents a superflow, if it is a local 
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Figure 20. Upper panel: Bifurcation of the trivial periodic state in a repulsive Gaussian comb 
for fi ~ 1. Lower panel: The unstable 27r-periodic wave function (dashed line) in comparison with 
the stable 47r-periodic wave functions (solid lines) in a Gaussian comb with A = 0.18 and fi = I. 



minimum of the Hamiltonian 

= / r (-^^ + Vix)^ ^ + - MlV'i'dx. (5.4) 

Otherwise a perturbation can lower the energy of the solution and the system becomes 
energetically unstable and superfluidity is lost, which is also called Landau-instability. 
Inserting the ansatz (j5.3p into the Hamiltonian and neglecting higher order terms reveals 
that the Bloch state is energetically stable if the operator 



with 



1 / d 



2 



^((?) = -2 l^^ + i^j +V{x) + 2g\xlj\'' - ^i. (5.6) 

is positive definite for all —1/2 < g < 1/2 [71]. 

Similarly one can deduce the dynamical stability by assuming a time-dependent per- 
turbation 

V'(x, t) = e'^^-'^' [uk{x) + 5uk{x, t)] . (5.7) 

The pertubation 6Ufi{x,t) grows exponentially in time if the operator azM{q) with cjz = 
diag(l, —1) has a complex eigenvalue for some q. Hence the Bloch state u^ix) is dynam- 
ically stable only if the spectrum of azM{q) is purely real for all —1/2 < q < 1/2. Note 
that dynamical instability always implies energetical instability, as shown in [71]. 

Here we only consider strictly periodic states, which are Bloch states with n = 0. 
The eigenvalues C{q) of azM{q) are calculated for the periodic states discussed above in 
dependence of the potential strength A. In Fig. [21] the maximum imaginary parts of the 
eigenvalues C,{q) are plotted, where the maximum is taken over all q £ [—1/2, 1/2]. One 
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Figure 21. Eigenvalues of azM{q), indicating the dynamical (in)stability, for the anti-symmetric 
(left) and symmetric (right) periodic solutions in dependence of the potential strength A. Shown is 
the maximum modulus of the imaginary parts of the eigenvalues of azM{q), where the maximum 
is taken over all eigenvalues and all q in the first Brillouin zone q G (—0.5, 0.5). 



observes that the spectrum of azM{q) is never purely real for the antisymmetric periodic 
states (except for A = 0). Hence the these states are always dynamically and energetically 
unstable. The symmetric periodic state is dynamically unstable for attractive potentials 
and not too strong repulsive potentials. It is stabilized if the potential strength A exceeds 
a critical value A^ > 0. However it still is energetically unstable for all values of A. Similar 
results are obtained for the delta-comb discussed in the previous sections. 

The dynamical (in) stability is illustrated in Fig. [22l where the time evolution of the 
squared modulus \ip{x,t)\'^ of a wavefunction is plotted. The initial state ■ijj{x,t = 0) is 
chosen as a symmetric periodic solution for A = 1 resp. A = 2 and = 1 plus a small 
random perturbation. According to the results from perturbation theory illustrated in 
Fig- l2H the system is dynamically unstable for A = 1 and dynamically stable for A = 2. 
This is well confirmed by the results of the wavepacket propagation. For A = 1 the onset of 
dynamical instability is clearly visible at t ~ 150, whereas no instability can be observed 
for A = 2. 

Note that the temporal stability of the states is not affected at all by the spatial stability 
properties. The onset of dynamical instability is solely determined by the spectrum of 
azM{q) and the strength of the initial random perturbation. On the contrary one finds 
that the parameter values for which a periodic state is spatially or dynamically stable 
are quite different. The lowest symmetric state is spatially stable only for Ac < A < 
(cf. Fig. [T2]) . whereas it is dynamically stable only for A > A^ > 0. (cf. Fig. [2T|) . One can 
even increase the potential strength A slowly through the critical value Ac for a bifurcation 
without affecting the onset of dynamical instability. 



6 Summary 



Summing up, the stationary solutions of the nonlinear Schrodinger equation for a delta- 
comb potential are analyzed. This somewhat artificial model system is of considerable 
interest since it allows analytical solutions in terms of Jacobi elliptic functions. 



28 



D Witthaut, K Rapedius and H J Korsch 




Figure 22. Time evolution of the squared modulus \ilj{x,t)\'^ of a wave function in a Gaussian 
comb of strength A = 1 (left) resp. A = 2 (right) The initial state i^{x,t = 0) was chosen as the 
symmetric periodic solution for = 1 plus a small random perturbation. 



The analogy to a classical nonlinear oscillator problem, the kicked nonlinear Hill equa- 
tion is pointed out and the dynamics in phase space is discussed. It is shown that a repul- 
sive nonlinearity typically leads to spatially chaotic dynamics and finally a divergence of 
the wave function, whereas the wave function always remains bounded for an attractive 
nonlinearity. Periodic solutions, which are fixed points of the phase space dynamics, are 
analyzed in detail. The fixed points are spatially stable only for certain parameter values, 
for which a stability map is calculated. The stability is lost in a Hamiltonian bifurcation 
scenario, giving rise to new period doubled fixed points. 

Nonlinear Bloch bands are calculated for the delta-comb showing the celebrated loop 
structures. It is shown that the Bloch states with zero quasi momentum arc just the 
periodic solutions mentioned above, which are given analytically in terms of Jacobi elliptic 
functions. The critical nonlinearity parameter for the emergence of loops is calculated 
analytically. Furthermore the emergence of 47r-periodic Bloch bands is discussed and 
linked to the period doubling bifurcation of periodic solutions mentioned above. 

Finally an extension to a more realistic potential, a periodic arrangement of narrow 
Gaussian potential barriers, shows that the properties of the stationary solutions remain 
qualitatively the same as for the delta-comb. In addition the temporal stability of the 
periodic solutions is analyzed. Antisymmetric periodic solutions are always dynamically 
unstable, whereas the symmetric periodic solutions are stabilized if the potential strength 
exceeds a critical value A > A,^ > 0. 

The properties of stationary states of the nonlinear Schrodinger equation are of funda- 
mental interest for the study of BECs in optical lattices. For example the emergence of 
loops leads to dynamical instability of BECs in tilted optical lattices because an adiabatic 
evolution is no longer possible here [25,39]. But not all solutions in periodic potential 
must take the form of Bloch states, as it was shown in the present paper. In contrast, 
most solutions of the nonlinear Schrodinger equation with a repulsive potential will be 
spatially chaotic or even divergent. This deserves further studies. 
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